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Abstract 

We consider a model consisting of a self-avoiding polygon occupying a variable density 
of the sites of a square lattice. A fixed energy is associated with each 90°-bend of 
the polygon. We use a grand canonical ensemble, introducing parameters \i and (3 to 
control average density and average (total) energy of the polygon, and show by Monte 
Carlo simulation that the model has a first order, nematic phase transition across a 
curve in the (3-fi plane. 
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1. Introduction 

In polymer physics self-avoiding walks have been used for many years to model 
long chain molecules [1,2]. One of the oldest such models, due to Flory [3], consists of a 
single self-avoiding random polygon occupying all the sites of a square lattice, with the 
randomness controlled by the total energy associated with 90°-bends in the polygon. 
Thinking of the polygon as made of many flexibly connected monomers, the model is a 
canonical ensemble with the temperature behavior analyzed only at the optimally high 
particle density of 1. At high temperature the Flory model behaves like a disordered 
fluid, while at low temperature the model displays long-range nematic order. (Flory 
was modeling the "melting" of a nematically ordered polymer [3].) Although it is 
generally accepted that there is a true phase transition between these regimes in the 
model, there has been a dispute over the character of the transition [4], which Flory 
had originally predicted to be first order. For the Flory model the controversy seems to 
have been settled by a recent paper [5] of Jacobsen and Kondev which shows that the 
transition is second order. (There is a useful summary of the history of the controversy 
in [5].) They also suggest, however, that the transition may become first order at 
lower density if vacancies are permitted (see Section VII C in [5]). In our paper we 
concentrate on this generalization of the Flory model which allows for vacant lattice 
sites, so that our random polygon is now controlled by two parameters, temperature 
and chemical potential. We show by Monte Carlo simulation that the model has a first 
order phase transition. We also use order parameters to illustrate the nematic nature 
of the transition. 

We note some other variations and applications of the Flory model. One variation 
uses an ensemble of polygons rather than one long polygon. Although the connection 
between the two models is unclear, it is of interest that there have been conflicting 
results on the character of the transition in this case also [6,7]. In a different direc- 
tion we note that variants of the Flory model have been applied to nonequilibrium 
materials, specifically to crumpled or confined sheets (in 3 dimensions) and wires (in 
2 or 3 dimensions). For instance in [8] we modeled a progressively confined wire in 2 
dimensions by a version of the model used here, but on a triangular lattice, showing 
a similar transition. A field theoretic (continuum) model of a variably confined wire 
by Boue and Katsev [9], however, found a second order transition, and it remains 
to understand the origins of this difference between their continuum model and our 
lattice models. 

2. The Model, and Results 

Consider the set W of all self-avoiding polygons (i.e., closed self-avoiding loops) 
on the square lattice with periodic boundary conditions, L = (Z/vZ) 2 , with v a fixed 
positive integer. The energy E(w) of polygon w in W is defined as the number of 
right angles in w, and the length N(w) of w is the number of unit line segments in w. 
Given an inverse temperature (3 = 1/T and a chemical potential //, the free energy of 
the model is j3E — (3(xN — S(E, N) where S(E, N) is the entropy, that is, the natural 
logarithm of the volume in phase space of self-avoiding polygons at fixed E and N. 
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Figure 1. a) The graph of average density vs. fx at j3 = 1.5, for system volumes 
V = 40 2 through V = 120 2 . b) Width of 95% confidence intervals for average density, 
for V = 120 2 . 

As usual in a grand canonical ensemble this is optimized by the probability mea- 
sure mp jfJ , defined on the subsets of W by 
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-/3(E(w)- fJ ,N(w)) 



1) 



for w G W, where Z$ ^ is the appropriate normalization. In this notation we have 
suppressed the dependence of mp^ and Zp tti on the system volume V = v 2 . 

To simulate the model we fix either (3 or fi, and then slowly increase the other 
parameter, starting from well into the disordered regime. The basic Monte Carlo step 
is as follows (see pgs. 41-44 in [10]). Given a polygon w(t) at step t in the simulation, 
we introduce, with probability pi, a trial configuration w(t)' which changes the length 
of w(t) by ii. If w(t)' is not self-avoiding then we take w(t + 1) = w(t); otherwise 
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we set w(t + 1) = w(t)' with probability q = min(Q, 1), and w(t + 1) 
probability 1 — q, where 



Q- 

Here pi = p 2 = 2/5, p 3 = 1/5, l\ 
a) 



f3[^i+E{w{t))-E{w{t)')} 



w(t) with 
2) 




Figure 2. a) The graph of average energy per volume vs. (3 at \x = 0.15, for system 
volumes V = 40 2 through V = 120 2 . b) Width of 95% confidence intervals for average 
energy per volume, for V = 120 2 . 

To determine whether the simulation for each pair (/?, fj.) has sufficiently many 
Monte Carlo steps, we compute a "mixing time" as the smallest t such that the stan- 
dard autocorrelation function 
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( _ 7s 2 y^^measjwi) - A) • (meas(w i+t ) - A) 
K n t)a , =1 
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falls below zero. Here meas represents any of our various measurements, described 
below, A and a 2 are the sample average and variance (respectively) of meas over the 
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simulation of (/3,fi), and Wi is the zth configuration in the simulation of /i), with 
n total steps. We found that our simulations of each (/3, fj.) were, in the worst cases, 
at least 5 mixing times long (on average), and we therefore believe our Monte Carlo 
runs are reasonably close to sampling the distributions mp^. We repeated each of our 
simulations 100 times, and obtained 95%- confidence intervals for meas from Student's 
t-distribution with 99 degrees of freedom on the average values of meas over each 
simulation. (Measurements related to specific heat were calculated differently, and 



are discussed below.) A single simulation of the largest system (V 
8 x 10 11 basic Monte Carlo steps. 



120 ) contains 




Figure 3. a) The graph of specific heat vs. /? at ji = 0.15, for system volumes 

V = 40 2 through V = 120 2 . b) Width of 95% confidence intervals for specific heat, for 

V = 120 2 . 



We measure average energy per volume {E)p^/V , average density (<^)/3, M , as well 
as order parameters corr and lay, which were introduced in [8] and are defined as 
follows. Given a polygon w, corr(w) is the proportion of edges in w which have the 
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same orientation (horizontal or vertical) as a randomly chosen edge in w. Given w, 
lay(w) is the normalized volume u 2 /V of the largest square sublattice L' = (Z/uZ) 2 C 
L such that the orientation of w (horizontal or vertical) at the origin agrees with the 
orientation of w at 80% or more of the sites in V . (We choose an orientation for the 
polygon w so that each lattice site has a unique horizontal or vertical orientation.) 
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Figure 4. a) The graph of specific heat vs. (3 and /U, for V = 100 2 . b) Estimation of 
/5 values which maximize specific heat at fixed \i. 

We compute specific heat (l/V)d(E) p^/dT from fluctuations, that is, 

T 2 ^^ = (E)pAvN - E)p,» - (E^N - E 2 )p^. 4) 
To compute the values of from equation (4), we took averages of the relevant 
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measurements meas over 100 independent simulations. Then for 95%-confidence in- 
tervals we repeated this process 4 times, and used the Student's t-distribution with 3 
degrees of freedom. We checked that the resulting curve agreed with the numerical 
derivative of energy. 

In contrast with [8] we simulate well into the ordered regime and find direct 
evidence of a first order phase transition. In particular the trends with increasing 
system size in the curves of Figs. 1 and 2 strongly suggest that the average density 
and average energy per volume both develop jump discontinuities at the transition, in 
the infinite volume limit. In confirmation, Fig. 3 shows the specific heat developing a 
delta function singularity at the transition. 



a) b) 
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Figure 5. a) The graph of corr vs. f3 at fi = 0.15, for volumes V = 40 2 to V = 120 2 . 
b) The graph of lay vs. /3 at \i = 0.15, for volumes V = 40 2 to V = 120 2 . 

We plot the specific heat surface as a function of f3 and /a, as well as the (/3, jj)- 
coordinates of the maximum of specific heat, at various 0.125 < \i < 0.175, 1.2 < (3 < 
1.5 in Fig. 4. The latter gives an indication of the transition curve; note that as \x 
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increases, the temperature at which the transition occurs increases. 

As evidence of an nematic transition, the measurements corr and lay (see Figs. 5 
and 6) exhibit a jump discontinuity at the transition from their disorder values of 1/2 
and zero, respectively. 
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Figure 6. a) The graph of corr vs. fx at f3 = 1.5, for volumes V = 40 2 to V = 120 2 . 
b) The graph of lay vs. /x at f3 = 1.5, for volumes V = 40 2 to V = 120 2 . 

3. Nonequilibrium 

We mentioned above that with variable density added, a version of the Flory 
model on a triangular lattice has been used [8] to model wires progressively confined 
in 2 dimensions. (See [11] for a related approach, and further references.) We add 
here a note on the modeling of such nonequilibrium materials. Since the polygons are 
self-avoiding, these lattice models introduce a unit length scale for the width of the 
wire. Now if it requires energy E to bend one wire of unit thickness to a given radius of 
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curvature it would require mE to bend a loose bundle of m parallel wires, but because 
of the interconnectedness it would require more than mE to bend a single wire of 
thickness m. Therefore the bending energy is highly nonlinear in the thickness of the 
wire, growing faster than the square of the thickness, and so energy is an independent 
parameter in our modeling. 




Figure 7. Configurations with V = 100 2 in equilibrium at: 

a) = (1.5,0.03) and density <f> = 0.392 

b) (/9,/i) = (1.5,0.13) and density <f> = 0.655, 

c) n) = (2.5, 0.01) and density <j> = 0.146, 

d) (13, n) = (2.5,0.025) and density </> = 0.392. 

Fig. 7 shows typical polygons with low and high volume fraction, but at two very 
different temperatures. This could be useful in estimating whether experimentally 
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confined materials are in an ordered or disordered regime. An alternative method 
would be to compute order parameters such as lay or corr. 

4. Conclusion 

We have introduced a version of the Flory model that allows for a positive fraction 
of vacancies, and shown by Monte Carlo simulation that the model has a first order, 
nematic phase transition. In terms of the (inverse) temperature (3 and chemical poten- 
tial \x of our grand canonical ensemble, we find that the transition lies approximately 
on the curve shown in Fig. 4 b). 

We conclude by contrasting the results in this paper with those in [8] , which uses 
a very similar model, but on a triangular lattice and therefore with greater complexity 
of the interaction energies. A grand canonical ensemble was also used in [8], but 
simulations were confined to a single isotherm. 

The evidence in [8] of a nematic transition is based on two order parameters, lay 
and corr, used also in this paper. The arguments in [8] are based heavily on trends in 
these order parameters as the system size is increased, namely that at low chemical 
potential fx both order parameters decrease monotonically toward their disordered 
values, while at higher fi both order parameters increase monotonically away from 
their disordered values. This strongly suggests a nematic transition, which would be 
expected to be first order. The simulations in [8] became unreliable moving into the 
ordered regime, so in particular no direct evidence was given of a discontinuity of a 
first derivative of the free energy, namely volume fraction or energy density. 

For the present paper, which also reports a nematic transition in a model similar 
to [8] but on a square lattice, as in the original Flory model [4,5], we were able to make 
reliable simulations well into the ordered regime, and simulated on a grid of values of \x 
and (3. The main improvement over [8] is that now we are able to show discontinuities 
in volume fraction and energy density, the usual hallmarks of a first order transition, 
as well as discontinuities in corr and lay, which clarify the nematic nature of the 
transition. We feel this is strong, direct evidence supporting the suggestion in [5] that 
the second order transition which they find at density 1 in the Flory model becomes 
first order at lower density when vacancies are included in the model. It still remains 
to reconcile this behavior with that found in [9] in their continuum model of confined 
loops. 
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